Numerical study of the coupled time-dependent Gross-Pitaevskii equation: 
Application to Bose-Einstein condensation 



Sadhan K. Adhikari 

Instituto de Fisica Teorica, Universidade Estadual Pauhsta, 01.405-900 Sao Paulo, Sao Paulo, Brazil 

(February 1, 2008) 



O 
O 

(N 



(N 



«4-H 

o 



o 



> 

in 
m 
o 



I 

O 

o 



X 



We present a numerical study of the coupled time-dependent Gross-Pitaevskii equation, which de- 
scribes the Bose-Einstein condensate of several types of trapped bosons at ultralow temperature 
with both attractive and repulsive interatomic interactions. The same approach is used to study 
both stationary and time-evolution problems. We consider up to four types of atoms in the study 
of stationary problems. We consider the time-evolution problems where the frequencies of the traps 
or the atomic scattering lengths are suddenly changed in a stable preformed condensate. We also 
study the effect of periodically varying these frequencies or scattering lengths on a preformed con- 
densate. These changes introduce oscillations in the condensate which are studied in detail. Good 
convergence is obtained in all cases studied. 
PACS Number(s): 02.70.Rw, 02.60.Lj, 03.75.Fi 



I. INTRODUCTION 

The experimental detection of Bose-Einstein con- 
densation (BEC) at ultralow temperature in dilute 
bosonic atoms (alkali and hydrogen atoms) employing 
magnetic traps have spurred intense theoretical activities 
on various aspects of the condensate . Many proper- 
ties of the condensate are usually described by the nonlin- 
ear mean-field Gross-Pitaevskii (GP) equation ^]^. The 
GP equation in both time-dependent and independent 
forms is formally similar to the Schrodinger equation in- 
cluding a nonlinear term [||J^:^- 

More recently, there has been experimental realization 
of BEC involving two types of atoms [^-|l2|. In the 
actual experiment ^^Rb atoms formed in the F — 1, 
m = — 1 and F = 2, m = 1 states by the use of a 
laser served as two different species of atoms, where F 
and m are the total angular momentum and its projec- 
tion In another experiment a coupled BEC was 
formed with the ^"^Rb atoms in the F — 1, m = — 1 and 
F = 2, m = 2 states. ||7|,p^. Experimentally, it is possi- 
ble to use the same magnetic trap to confine the atoms 
in two quantum states, which makes this study easier 
technically compared to the formation of a BEC with 
two different types of atoms requiring two different trap- 
pingmcchanisms. It has also been found in these studies 
|jll| , [l2| that the ®^Rb atoms have a repulsive interaction 
in all three states considered above. These experiments 
initiated theoretical activities in multicomponcnt BEC 
described by the coupled GP equation . 

A numerical study of the time-dependent coupled GP 
equation is interesting as this can provide solution to 
many time-evolution problems involving more than one 
types of atoms forming a BEC. The solution of the cou- 
pled nonlinear GP equation is nontrivial [ p^ and here wc 
undertake the challenging task of the numerical study of 
these time-evolution problems. 

In a multi-component BEC the main feature is the cou- 



pling between different types of condensates, which can 
lead to new effects associated with this novel and surely 
much richer situation inexistent in a single-component 
BEC. We list some interesting possibilities below which 
we shall investigate numerically in this study. It is possi- 
ble to have a distinct trapping frequency for each of the 
components, one of which can be suddenly altered exper- 
imentally. Of course, this change would affect the compo- 
nent of the BEC directly trapped by this field. However, 
it is more interesting to study how this sudden change 
affects the other component of BEC not directly trapped 
by this field. Also it is possible to suddenly vary jl7| the 
atomic scattering length of one of the species and study 
its effect on the other component. The above variation 
of one of the trap frequencies or scattering lengths can 
be carried through in a periodic fashion and its effect on 
the other component can be studied. The studies men- 
tioned above are peculiar to a coupled BEC and are of 
interest as it is now possible to vary both the trap fre- 
quencies 1 16 as well as the scattering lengths [0,|l 17 
both abruptly or in a periodic fashion. These effects do 
not have any analogue in the uncoupled BEC and we 
motivate the present investigation with special empha- 
sis on these effects in this study. We investigate these 
time-evolution problems using a set of two coupled GP 
equations in the purely repulsive case. In addition we 
study the stationary solution to the coupled GP equa- 
tion describing a multi-component condensate, where we 
consider up to four components. 

We solve the coupled BEC problem using the time- 
dependent coupled GP equation in cases of attractive and 
repulsive atomic interactions by discretization with the 
Crank-Nicholson-type rule complimented by the known 
boundary conditions at origin and infinity [ p^ . This pro- 
cedure leads to good convergence for both the stationary 
and time-evolution problems. 

First we consider stationary coupled condensates un- 
der the action of trap potentials. Stable and converged 
numerical results are obtained for up to four coupled 
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equations in the repulsive case and two in the attrac- 
tive case. The time-dependent GP equation is directly 
solved to obtain the full time-dependent solution in the 
case of stationary problems, from which a trivial time 
dependent phase factor is separated and the stationary 
solution obtained as in the uncoupled case § . 

We also study the numerical stability of the calcula- 
tional scheme which is more difficult to obtain when the 
nonlinearity is large. For this purpose we only consider 
repulsive interatomic interactions in the time-evolution 
problems where condensates with large nonlinearity can 
be formed. In the case of attractive interaction a large 
nonlinearity leads to the collapse of the condensate Q . 

In Sec. II we present the coupled time-dependent GP 
equation which we use. In Sec. Ill we describe the nu- 
merical method in some detail. In Sec. IV we report the 
numerical results for the stationary case and in Sec. V we 
report a study of three types of time-evolution problems. 
Finally, in Sec. V we give a summary of our investigation. 



II. NONLINEAR COUPLED GROSS-PITAEVSKII 
EQUATION 

The GP equation for a coupled trapped Bose- 
Einstein condensate at zero temperature is written as 



1 ^ 

— V2 -I- -cmwV^ +y2g,iNi\-^i{r, r)p 



— ih 



dr 



1=1 

^',(r,T) =0, 



(2.1) 



where ^i(r, t) at position r and time r is the wave func- 
tion for the component i of the condensate, m is the 
mass of a single bosonic atom, Ni the number of con- 
densed atoms of type I, M the number of types of atoms, 
Cjmw^r^/2 the attractive harmonic-oscillator trap poten- 
tial, uj the oscillator frequency. The parameter Ci has 
been introduced to independently modify the frequency 
of the harmonic oscillator trap for each type of atoms. 
Here gu — 4,7111^ an /m is the coupling constant for elastic 
interaction between atoms of types i and I, where an is 
the corresponding scattering length. The masses of dif- 
ferent types of atoms are taken to be equal, as this is 
necessary while considering the coupled BEC formed of 
different spin states of ^^Rb — one of the most important 
realization to date. In this work we shall not allow the 
transition of one type of atoms of the BEC to the other 
and take the number of atoms of each component to be 
constant as in the experiment of Ref . . 

Here we shall be interested in the spherically symmet- 
ric solution 5'i(r, r) = ipii^r, r) to Eq. ( ^.l| ), which can be 
written as 
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dr 



V'i(r,r) = 0, 



(2.2) 



The above limitation to the spherically symmetric solu- 
tion (in the zero angular momentum state) reduces the 
coupled GP equation to a one-dimensional coupled par- 
tial differential equation. 

As in Refs. it is convenient to use dimcnsionless 
variables defined by x = ^/2r/a]^Q , and t — tlo, where 
aj^Q = yJn/{mLu), and 0j(a;,t) = xipi{r,T){-\/2iia^^Y/'^ . 
In terms of these variables Eq. becomes 
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dt 



1=1 

4>i{x, t) = 0, 



(2.3) 



where nn = 2^/2Niaii / a-^^ is the reduced number of par- 
ticles and this number could be negative (positive) when 
the corresponding scattering length is negative (posi- 
tive) representing an attractive (repulsive) interatomic 
interaction. The normalization of the wave function is 
/o°° l^'iix J d,x — 1 and its root-mcan-squarc (rms) ra- 



dius Xp*j^g is given by 



(i) 

^rms 



-,1/2 



x'^\(j)i{x, t)\'^dx 



(2.4) 



III. NUMERICAL METHOD 



To solve Eq. (2.3) numerically one needs the proper 
boundary conditions at x — > and cxi. For a confined 
condensate, for a sufficiently large x, (f>i{x,t) must van- 
ish asymptotically. Hence the nonlinear term propor- 
tional to \(j)i{x,t)\^ can eventually be neglected in the 
GP equation for large x. Consequently the the asymp- 
totic form of the physically acceptable solution is given 
by lima;^oo |0i(a^, i)| ~ exp(— a;2/4). Next we consider 
Eq. ( |2.3| ) as x — > 0. The nonlinear term approaches 
a constant in this limit because of the regularity of the 
wave function at x = 0. Then one has the condition 
10,(0, i)| =0. 



A convenient way to solve Eq. (2^) numerically is to 
discretize it in both space and time and reduce it to a 
set of algebraic equations which could then be solved by 
using the known asymptotic boundary conditions. The 
procedure is simil ar t o that in the uncoupled case 
We discretize Eq. (2.3) by using a space step h and time 
step A with a finite difference scheme using the unknown 
{(f>i)p which will be approximation of the exact solution 

4'i{^P^ tk) where Xp = ph and tk = kA. The time deriva- 
tive in Eq. (2.3) involves the wave function at times tk 
and ifc+i. As in the uncoupled case we expr ess the wave 
functions and their derivatives in Eq. (^.3[) in terms of 
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the average over times tk and tfc+i and the resultant 
scheme leads to accurate results and good convergence. 
In practice we use the following Crank-Nicholson-type 
scheme |l5| to discretize the partial differential equation 
(PI) 



h)l 



A 



1 

2^ 



0^11 -2( 



\k+l 



1 
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Ci Xp 
\k+l 



)p-i 

pi 



1=1 



(3.1) 



Considering that the wave function components (pi are 
known at time tfc, Eq. ( |3.1| ) is an equation in the un- 
knowns - ('/>i)p+i: (0Op"^^ and {4>i)lt\- In a lattice 
of N points Eq. (3.1) represents a tridiagonal set for 
p = 2, 3, (iV — 1) for a specific component This set 
has a unique solution if the wave functions at the two end 
points {4'i)\'^^ and (00^^^ are known ||l^. In the present 
problem these values at the end points are provided by 
the known asymptotic conditions. 

To find the ground state of the condensate we start 
with the analytically known properly normalized wave 
functions of the uncoupled harmonic oscillator problems 
described by Eq. ( |2.3| ) with nu = 0. We then repeat- 
edly propagate these solutions in time using the Crank- 
Nicholson- type algorithm (3.1). Starting with nu = 0, 
at each time step we increase or decrease the nonlinear 
parameter nu by an amount Ai. This procedure is con- 
tinued until the desired final value of nu is reached. The 
resulting solution is the ground state of the condensate 
corresponding to the specific nonlinearity. 

The time-dependent approach is the most suitable for 
solving time-evolution problems. In the present study 
we only consider evolutions problem starting from a sta- 
ble condensate at t = 0. In these cases the stationary 
problem is solved first and the wave function so obtained 
serves as the starting wave function for the time-evolution 
problem. 



IV. RESULTS FOR THE STATIONARY 
PROBLEM 

First we consider the stationary ground-state solution 
of Eq. (2_^) in cases of both attractive and repulsive in- 
teractions using two and four coupled equations. The 
numerical integration was performed up to xniax = 15 
with h = 0.0001 using time step A — 0.05 and the pa- 
rameter Ai = 0.01. After some experimentation we find 
that good convergence is obtained with parameters A 
and Ai near these values. The convergence is fast for 
small nonlinearity. The final convergence of the scheme 
breaks down if nonlinearity is too large. In practice these 



difficulties start for nu > 20 for the ground state for a re- 
pulsive interaction in a computational analysis in double 
precision. For an attractive interaction the coupled GP 
equation does not sustain a large nonlinearity \nii\ and 
leads to collapse. Except for values of nonlinearity near 
collapse, the GP equation in the attractive case leads to 
good convergence. 



A. Repulsive Atomic Interaction 

In most of the experimental realization of BEG in 
trapped atoms, the interatomic interaction is repulsive 
and we consider this case first. We consider the simple 
case of two coupled GP equations in the case of repulsive 
interaction with (a) nu = 7122 = 10, ni2 = 71.21 = 5, 
ci = 1 and C2 = 0.25, (b) nu = ^22 = 5, ni2 = n2i = 
100, ci = 1 and C2 — 0.25, and (c) nu = n22 = 10, 
ni2 — n2i = 5, and ci = C2 = 1. In this case all 
interactions are repulsive corresponding to the positive 
sign of all nu = 2\/2Niau / cl-^q- Although these param- 
eters are in dimensionless units, it is easy to associate 
them to an actual physical problem of experimental in- 
terest. For example, for the mixture of |F = l,m = — 1) 
and |F = 2,rn = 1) ^^Rb states, the ratio of scatter- 
ing lengths a|i,_i)/a|2,i) = 1.062 ||l|]. If we label the 
|1, —1) state by 1 and the |2, 1) state by 2, and consider 
aii/aj^Q — 022/0^0 — 0.002, then nu — ^22 = 10 cor- 
responds to A^i ~ A^2 — 1770. This estimate gives an 
idea of the actual experimental condition that the present 
set of parameters simulate. The three models considered 
above can simulate actual experimental situations com- 
posed of two states of ^^Rb. The different values of ni2 
and n2i considered above corresponds to different possi- 
ble unknown repulsive interactions among the two species 
of condensates. It is realized from the coupled GP equa- 
tion that in case of model (c) (pi — (j)2- We show results 
for the two components of the wave function for two sets 
of values of Ai: 0.01 (full line) and 0.015 (dashed line) 
in Fig. 1. The difference between the two sets of results 
increases as the nonlinearity of the GP equation given 
by nu increases, e.g., for the case (b) above compared to 
(a). The difference reduces to zero as the nonlinearity 
decreases. 

Next we consider the more complicated case of four 
coupled GP equations with repulsive atomic interaction. 
This is a purely theoretical case with no experimental 
analogue, as all experiments to date are limited to two 
coupled condensates only. In this case the numerical 
method works in the same fashion and good convergence 
is obtained with moderate values of nonlinearity. Again 
we consider the case of repulsive interaction between all 
possible pairs. In this case in the four-component model 
we take nu — 4, n22 = 5, n^s = 6, n^^ — 8, and 
nu = 2, z 7^ ci = 4, C2 = 1, C3 = 4, and C4 = 0.25. 
The solution for the wave-function components obtained 
with A = 0.05 and Ai = 0.01 (full line), and Ai = 0.015 
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(dashed line) are shown in Fig. 2. The maximum dif- 
ference between the two calculations is about 6% for the 
largest component (0i) near a; = 0, although the average 
difference is much smaller as can be seen in Fig. 2. 

B. Attractive Atomic Interaction 

The case of attractive atomic interaction demands spe- 
cial attention as one can have the phenomenon of col- 
lapse in this case. We consider the case of two cou- 
pled GP equations with attractive interactions between 
like atoms ii, i = 1,2, and with repulsive interactions 
between unlike atoms ij, i ^ j. In this case some of 
the atomic interactions are attractive in nature. Con- 
sequently, with large (attractive) nonlinearity the sys- 
tem may undergo collapse and for stable stationary solu- 
tion of the GP equation the nonlinearity should be main- 
tained small. We consider the following three cases: (a) 
nil = -1, "22 = -1.5, ni2 = "21 = 4, ci = 4, and 
C2 = 1; (b) nil = -1, "22 = -1.5, ni2 = "21 = 80, 
ci = 4, and C2 = 1; (c) nn = -1, ^22 = -1.5, 
"12 — "21 = 100, ci = 4, and C2 = 1. The only ob- 
served case of BEC with attractive interaction is the case 
of "^Li with a/aj^Q ~ -0.0005 0. If we label this state 
by 1 then nn — —1 corresponds to the actual particle 
number A^i = 700. In an uncoupled condensate of ^Li 
the BEC collapses for more than 1400 atoms. Although 
there has been no experimental realization of coupled 
BEC in the case of attractive interaction, the parameters 
cited above may simulate the BEC of ground-state atoms 
of ^Li coupled to one of its excited states, where the 
atomic interaction is also attractive. If we assume that 
the excited-state atoms have the same value of a/aj^o ^® 
in the ground state then 7122 = — 1.5 corresponds to the 
number of atoms N2 = 1500 in the excited state where 
the excited state is labeled by the index 2. 

The wave-function components in this case are shown 
in Fig. 3. As for stable stationary solution the nonlin- 
earity in this case has to be smaller than in the purely 
repulsive case, numerically it is easier to obtain precise 
solution except for values of the nonlinearity close to (and 
beyond) collapse. The parameters above in cases (a), (b), 
and (c) are chosen to illustrate the collapse of the sys- 
tem arising from the divergence of the first component 
of the wave function ((/fi). The nonlinear couplings ni2 
and 7121 are increased as we move from case (a) to (b) 
and then to (c) . Other parameters of the model are kept 
fixed. Because of the attractive interactions in this case, 
the system moves towards collapse as we move from case 
(a) to (c) through (b). The central density correspond- 
ing to (/)i increases (eventually tends to infinity) and the 
rms radius of the system decreases (eventually tends to 
zero) with the increase of nonlinearity. This is clear from 
a comparison of Fig. 3 with Figs. 1 and 2. The very 
large central value of the wave function (pi (~ 35) and 
its small radial extension indicates a large central density 



and a small rms radius. 

Finally, we consider the case of two interacting systems 
with all interactions repulsive. In this case the system is 
more vulnerable to collapse if the nonlinearity is large. 
We consider the following three sets of parameters in this 
case for which we show the solution in Fig. 4: (a) nu = 
-1, 7122 = -1, "12 = "21 = -0.4, ci = 4, and C2 = 0.25; 
(b) 7iii = -1, 77,22 = -1, "12 = "21 -0.5, ci = 4, and 
C2 = 0.25; (c) nil = -1, "22 -1, "12 = "21 = -0.55, 
ci = 4, and C2 = 0.25. These parameters simulate the 
possible coupled BEC composed of the attractive ground 
and excited states of ^Li where the interaction between 
a ground- and an excited-state atom is also taken to be 
attractive. It is possible to calculate the number of the 
two types of atoms as in the discussion related to Fig. 
3. Here the nonlinearity increases as we move from case 
(a) to (b) and then to (c), and consequently, the wave- 
function components become more and more localized 
with a large central density and small rms radii signaling 
the onset of collapse of the system. This is clear from 
Fig. 4. In case (c) the nonlinearity is the highest and 
one is closer to collapse. However, there is a difference 
between the two collapses shown in Figs. 3 and 4. In Fig. 
3 the route to collapse is manifested through a singular 
behavior of component (pi of the wave function; in Fig. 4 
both components <pi and (p2 exhibit the singular behavior. 
The spacial extension of the wave function components 
close to collapse in Figs. 3 and 4 is much smaller than the 
wave function components in the purely repulsive cases 
shown in Figs. 1 and 2. 

C. Estimate of Numerical Error 

It is appropriate to comment quantitatively on the nu- 
merical accuracy of the present method. If we iterate 
the final solution in time without changing the nonlin- 
earity, the numerical result keeps on oscillating with a 
small amplitude around the converged value. This oscil- 
lation gives a good estimate of the numerical error of the 
method. This error manifests in a different way in Fig. 
1, where we have varied Ai. The numerical solution of 
the time-dependent method is independent of the space 
step h provided that a typical value around h = 0.0001 
is employed as in the present study. No visible difference 
in the solution is found if h is increased by a factor of 
two. However, the above oscillations with time iteration 
are sensitive to the parameters A and Ai. The values 
of these parameters (A = 0.05 and Ai = 0.01) are cho- 
sen to minimize the oscillation of the results with time 
iterations. The oscillation increases if larger or smaller 
values of one or both of these parameters are employed 
and can really be large if an improper value of step A or 
Ai is chosen. This oscillation is quite similar to that in 
the uncoupled case studied in detail in Ref. ||]. 

From Figs. 1 and 2 we find that the error in the wave 
function component \(l)i{x,t)/x\ as a function of x is the 
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largest at a; = 0. In Table I we show the percentage er- 
ror of lim^^o \(l)i{x,0)/x\, defined by £ = 100[|(?iii(0, t)| — 
|</'i(0, 0)|]/|(^i(0, 0)1 at those iterations where this error 
is maximum. For illustration we consider the models 
(a) discussed in Figs. 1 and 4. At positive £ there is 
overshooting and at negative £ there is undershooting. 
Between a positive £ and a negative £ there is a zero 
of £ denoting zero error. We find that the wave func- 
tions oscillate with time around the stationary solution. 
The maximum reported error is about 6%. Consider- 
ing that we are dealing with coupled nonlinear equations 
these errors are well within the acceptable limit. The 
errors shown in Table I would also be the typical errors 
in time-evolution problems with same nonlinearity which 
we study in Sec. V. From Table I we see that the pe- 
riod of oscillation of the result varies from one case to 
another. However, the error increases as the nonlinearity 
increases or as the system approaches collapse in the case 
of attractive interaction. For example, in models (c) of 
Figs. 3 and 4, which are close to collapse corresponding 
to almost maximum permissible nonlinearity, the error 
increases quickly with time iteration and a large numer- 
ical error could be generated easily. 

V. RESULTS FOR THE TIME-EVOLUTION 
PROBLEM 

Now we consider three types of time-evolution prob- 
lems, some of which could possibly be studied experimen- 
tally. As the repulsive case leads to more stable configu- 
ration of the condensate, in this study we consider only 
this case in the process of time evolution. The attractive 
case of coupled BEC is also very interesting from a phys- 
ical point of view because of the occurrence of collapse. 
We have performed a study of the dynamics of collapse in 
coupled BEC using the present numerical method, which 
will be reported elsewhere. The two types of parameters 
that can be varied in the time-evolution study are the 
frequencies of the harmonic oscillator traps and the dif- 
ferent scattering lengths. Recently, it has been possible 
to vary the scattering length experimentally by varying 
an external field |pl] , |l2| . It is also possible to vary the 
trap frequency by varying the currents in the magnets 
responsible for confinement Jl6t . 

In the first type of problems we consider a sudden 
change of the harmonic oscillator frequencies or scatter- 
ing lengths at t = and study its effect on a preformed 
condensate. In the second type we study the effect of 
a periodic temporal variation of these frequencies on a 
preformed condensate. Finally, we study the effect of a 
periodic temporal variation of the scattering lengths on 
the preformed condensate. In all cases we take the pre- 
formed condensate as the one described by the model (a) 
of Fig. 1. We have commented before that the param- 
eters of this model can simulate the coupled BEC com- 
posed of the ground and an excited state of *^Rb, which 
gives a motivation for this choice. When we implement 



these time-dependent perturbations, the system starts to 
oscillate (grow and shrink) with time. The corresponding 
evolution can be studied best through the rms radii Q] 
which execute periodic oscillation with time. 

A. Sudden Change of Trap Frequency or Scattering 
Length 

By varying the external fields it is possible to vary the 
harmonic oscillator trap frequency of the confining trap 
as well as the atomic interactions (scattering lengths) 
]l7t . First we consider a sudden change of both the trap 
frequencies on the preformed coupled stationary BEC 
state corresponding to model (a) of Fig. 1. We set the 
reduced time T = t/0.05 = when we start the time evo- 
lution. As the time step A is 0.05, T is just the number 
of iterations. In this model we have two different trap 
frequencies for the two components given by ci = 1 and 
C2 = 0.25. In the first case at T = we suddenly inter- 
change the constants ci and C2, i.e., we set ci = 0.25 and 
C2 = 1 and study the time evolution. The evolution of 
the rms radii corresponding to 0i and 02 sue shown in 
Fig. 5 (a). Both rms radii execute oscillations. However, 
that corresponding to has a much larger amplitude. 
The periods of oscillation of the two radii are different. 

Next we consider a sudden change in one of the trap 
frequencies on the preformed condensate at T = corre- 
sponding to model (a) of Fig. 1 with ci = 1,C2 = 0.25. 
At T = 0, we set ci = C2 = 1, which corresponds to pa- 
rameters of model (c) of Fig. 1 . The evolution of the rms 
radii is shown in Fig. 5 (b). Although in the stationary 
configuration of model (c) of Fig. 1, = 02, this con- 
dition is never attained in this evolution problem. The 
system keeps on oscillating indefinitely with time. The 
oscillation shown in Figs. 5 (a) and (b) has nothing to 
do with the nonlinear or coupled nature of the problem. 
Similar oscillation also appears in an uncoupled linear 
oscillator when the trap frequency is suddenly changed. 
In the present coupled nonlinear problem both rms radii 
execute oscillations with time. However, when the am- 
plitude of oscillation of one of the components increases, 
that of the other decreases. This behavior denotes the 
transfer of kinetic energy from one component to the 
other. 

Now we study the effect of a sudden change of the scat- 
tering length(s) on the preformed condensate We 
consider the problem when the parameters of model (a) 
of Fig. 1 are suddenly changed to those of model (b) of 
Fig. 1 at T = 0. This is achieved by changing the non- 
linearities suddenly at T = from nn ~ n22 — 10, ni2 — 
rt2i = 5 to nil = — 5,rti2 = n2i — 100 with a vari- 
ation of the external field which controls the scattering 
length(s). In this case the oscillations of the system are 
shown in Fig. 5 (c) where we plot the time evolution of 
the rms radii of the two components. Both components 
of the condensate execute oscillations but with different 
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frequencies and amplitudes. One of the components exe- 
cute giant oscillations with large amplitude, whereas the 
amplitude of the other is much smaller. 

Finally, we consider the case when one of the trap- 
ping potentials is switched off at T = on the preformed 
condensate of Fig. 1 model (a) by setting C2 = 0. The 
oscillation in this case is shown in Fig. 5 (d) where we 
plot the two rms radii. In the absence of the trapping 
potential the second component of the condensate can 
no longer remain localized in space. However, it does not 
expand monotonically before evaporating. It starts to ex- 
ecute giant oscillation and eventually escapes to infinity. 
Similar oscillation was found in the case of an uncou- 
pled BEG when the trapping potential was removed . 
The first component essentially remains unchanged dur- 
ing the process under the action of the unchanged trap 
potential. The minor oscillation of the rms radii of the 
first component is due to the coupling to the expanding 
second component. 

B. Periodic Oscillation of Trap Frequency 



nil = 1 — 0.5 sin(7rT/20) for T > on the preformed con- 
densate of model (a) of Fig. 1. The resultant oscillation 
of the rms radii are shown in Fig. 7. This variation cor- 
responds to a variation of the atomic interaction among 
atomic states of the first type. Consequently, the rms 
radii of the first component of the BEC executes pro- 
nounced oscillation with moderate amplitude. There is 
no direct variation in the parameters of the second com- 
ponent. The second component of the condensate feels 
the effect of variation of rin through the coupling to the 
first component. Because of this secondary effect the sec- 
ond component also executes oscillation as can be seen 
from its rms radii in Fig. 7, albeit with a much smaller 
amplitude compared to the first component. 

We also considered a periodic variation of the scatter- 
ing length between one atom of each type (ai2 = 021) by 
setting ni2 = n.21 — 0.5 — 0.25 sin(7rT/20) on the same 
preformed condensate for T > and studied the resul- 
tant oscillation of the rms radii. However, no interesting 
behavior was observed and we do not show details of that 
oscillation. 



Instead of making a sudden change in the parameters 
of the model, next we introduce periodic oscillation in 
some of the parameters of the model for T > and study 
the consequence on the system. We introduce a peri- 
odic variation in the parameters Ci which are related to 
the harmonic oscillator trap frequencies. Experimentally, 
this variation is possible via a variation of the external 
fields which are controlled by currents. 

We again consider at T = the preformed condensate 
of the model (a) of Fig. 1. First we consider the variation 
ci = 1 — 0.5 sin(7rT/20) which corresponds to varying the 
frequency of the first trap. The resultant variation of the 
two rms radii are shown in Fig. 6. The first radius (full 
line) oscillates more rapidly with larger amplitude and 
frequency than the second radius (dashed line). This is 
reasonable as we are directly varying the first frequency 
in this case. The rms radius of the second wave function 
feels the effect through its coupling to the first compo- 
nent. We also varied both the papameters ci and C2 in a 
periodic fashion which corresponds to varying both the 
frequencies. In this case both the rms radii execute os- 
cillation. However, no interesting effect is observed and 
we do not show the details of that oscillation here. 



C. Periodic Oscillation of Scattering Length 

Now we study the effect of a periodic variation of the 
scattering length(s) of the system on a preformed con- 
densate. In our formulation this corresponds to a peri- 
odic variation of the parameters nu. This variation of 
the atomic interactions or the scattering lengths is now 
feasible experimentally We consider the periodic 

variation in one of the scattering lengths (an) by setting 



VI. SUMMARY 

In this paper we present a numerical study of the cou- 
pled time-dependent Gross-Pitaevskii equation for BEC 
in three space dimensions under the action of harmonic 
oscillator trap potentials with attractive and repulsive in- 
terparticle interactions between different types of atoms 

ill. 

The time-dependent coupled GP equation is solved 
by discretizing it using a Crank-Nicholson-type scheme 
|H,|5|. This method leads to good convergence for small 
nonlinearity. Numerical difficulty appears for large non- 
linearity (nu > 20). For medium nonlinearity, the ac- 
curacy of the method can be increased by reducing the 
space step h. 

The ground-state stationary wave functions are found 
to be sharply peaked near the origin for attractive inter- 
atomic interaction for larger nonlinearity (Fig. 4). For 
a repulsive interatomic interaction the wave function ex- 
tends over a larger region of space (Figs. 1 and 2). In 
the case of an attractive potential, the rms radii decrease 
with an increase of nonlinearity. There could be a col- 
lapse for attractive interaction when the nonlinear pa- 
rameters nu are increased as in the uncoupled case [^j. 
In the purely repulsive case we solved two and four cou- 
pled GP equations. In problems involving attraction we 
solved only the two coupled GP equations. 

In addition to the stationary problem we studied three 
types of evolution problems. A stable coupled condensate 
is considered at T = on which a time-dependent per- 
turbation is introduced. Two types of perturbations were 
considered on a two-component condensate with purely 
repulsive interactions. In the first type a sudden change 
in the parameters related to the frequencies of the trap 
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and the scattering lengths was introduced. In the sec- 
ond type a periodic variation of the different scattering 
lengths and the frequencies of the harmonic oscillator 
trap was introduced. In all cases the condensates ex- 
ecute periodic oscillation which is studied via the time 
evolution of the rms radii as in the uncoupled case 0. 
We conclude that the present time-dependent approach is 
very suitable for studying both the stationary and time- 
evolution problems of a coupled BEC. 

The work is supported in part by the CNPq and 
FAPESP of Brazil. 
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Figure Caption: 

1. Wave function components (pi (label 1) and (/)2 (la- 
bel 2) for two coupled GP equations with (a) rin = 
«22 = 10, ni2 = n2i = 5, ci = 1, C2 = 0.25; (b) 
nil = n22 = 5, ni2 = n2i = 100, ci = 1, C2 = 0.25; and 
(c) nil = «22 = 10, ni2 = n2i = 5, ci = 1, C2 = 1 calcu- 
lated with A = 0.05 and Ai = 0.01 (full line); A = 0.05 
and Ai = 0.015 (dashed line). In case (c) only the results 
for Ai = 0.01 arc shown. 

2. Wave function components (/>i (label 1), ^2 (label 
2), (^3 (label 3), and (^4 (label 4) for four coupled GP 
equations with riu = 4, n22 = 5, 7133 = 6, 7144 = 8, and 
nu = 2, i ^ I; ci = 4, C2 = 1, C3 = 4, and C4 = 0.25 
calculated with A = 0.05 and Ai = 0.01 (full line); A = 
0.05 and Ai = 0.015 (dashed line). 

3. Wave function components (j)i (label 1) and (j)2 (la- 
bel 2) for two coupled GP equations with (a) nu = — 1, 
n22 = —1.5, ni2 = n2i = 4, ci = 4, and C2 = 1 (dashed- 
dotted line); (b) nu = —1, n22 = —1.5, ni2 = n2i = 80, 
Ci = 4, and C2 = 1 (dashed line); (c) nu = — 1, 
«22 = -1.5, ni2 = n2i = 100, ci = 4 , and C2 = 1 
(full line) calculated with A = 0.05 and Ai = 0.01. 

4. Wave function components ipi (label 1) and (f>2 (la- 
bel 2) for two coupled GP equations with (a) nu = — 1, 
n22 = —1, 7^12 = ^21 = —0.4, ci = 4, and C2 = 0.25 
(dashed-dotted line); (b) nu = —1, 7122 = —1, ^12 = 
^^21 = —0.5, ci = 4, and C2 = 0.25 (dashed line); (c) 
nil = -1, "22 = -1, ni2 = n2i = -0.55, ci = 4, 
and C2 = 0.25 (full line) calculated with A = 0.05 and 
Ai = 0.01. 

5. The rms radii of the two components (f>i (full line) 
and (p2 (dashed line) of the wave function at different 
reduced times T = t/0.05 for the oscillating condensate 
when on the preformed condensate of model (a) of Fig. 
1 we suddenly inflict the following changes: (a) ci = 
0.25 and C2 = 1; (b) ci = C2 = 1; (c) nn = n22 = 5; 
^21 = ni2 = 100; (d) C2 = calculated with A = 0.05 
and Ai = 0.01. All other parameters are maintained 
unchanged. 

6. The rms radii of the two components ^1 (full line) 
and (1)2 (dashed line) of the wave function at different 
reduced times T = i/0.05 for the oscillating condensate 
when on the preformed condensate of model (a) of Fig. 
1 we suddenly inflict the following change: Ci = 1 — 
0.5sin(7rT/20) calculated with A = 0.05 and Ai = 0.01. 
All other parameters are maintained unchanged. 

7. The rms radii of the two components (f>i (full line) 
and (j)2 (dashed line) of the wave function at different 
reduced times T = i/0.05 for the oscillating condensate 
when on the preformed condensate of model (a) of Fig. 
1 we suddenly inflict the following change: nn = 1 — 
0.5sin(7rT/20) calculated with A = 0.05 and Ai = 0.01. 
All other parameters are maintained unchanged. 
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Table 1: Percentage error £ = 1OO[|0,(O, r)| - |(/),(0, 0)|]/|(/)j(0, 0)| of \<j).,{Q,T)\ [i = 1,2) at successive reduced times 
r(= i/0.05), where this error is maximum, calculated with A = 0.05 and Ai = 0.01. The cases considered correspond 
to model (a) of Figs. 1 and 4 
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